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Abstract 

The detection of insufficiently resolved or ill-conditioned integrand structures is 
critical for the reliability assessment of the quadrature rule outputs. We discuss a 
method of analysis of the profile of the integrand at the quadrature knots which 
allows inferences approaching the theoretical 100% rate of success, under error esti- 
mate sharpening. The proposed procedure is of the highest interest for the solution 
of parametric integrals arising in complex physical models. 
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1 Introduction 



A large number of physical models currently under study are characterized by 
two combined features. First, the observables are obtained as integrals which 
cannot be solved analytically. Second, the models describe physical systems 
involving one or more specific parameters the variation of which results in 
critical modification of the system behaviour. As a consequence, deep under- 
standing of the predictions of the models needs the exploration of the values 
of the observables over a large range of the variable parameters. 

1 Corresponding author; e-mail: adamg@theory.nipne.ro 

2 e-mail: adams@theory.nipne.ro 

3 e-mail: plakida@thsunl.jinr.ru 



Preprint submitted to Elsevier Preprint 



1 February 2008 



As usual, to solve the occurring parametric integrals, recourse is made to ex- 
isting library codes of automatic adaptive quadrature which may fail badly 
without providing any hint about such possibilities. We are directly aware 
of three such frustrating experiences. The first one concerns the two-band 
singlet-hole Hubbard model of cuprate superconductors [l]-[3], which involves 
integrals over ranges of the first Brillouin zone. The variation of the parameter 
of the model (the hole or electron doping in the high-T c superconductor) re- 
sults in substantial modification of the behavior of the involved functions over 
the Brillouin zone. The exploration of the predictions of the physical model 
with the doping is fundamental for the validation of the proposed mechanism 
as responsible for the superconducting pairing in cuprates. However, the relia- 
bility of the outputs was found to be exceedingly low to allow sound inferences 
based on the bare numerical outputs. A similar problem arises in the alterna- 
tive U(l) x SU(2) gauge theory model of underdoped cuprate superconductors 
[4]. The meaningful physical solution derived under simplifying assumptions in 
[5] could not be recovered from outputs generated by the available automatic 
adaptive quadrature codes. The numerical exploration of a model of nuclear 
fission [6] could not be achieved by means of library quadrature codes either. 

These circumstances come from the fact that the existing algorithms for the 
numerical integration of real valued functions (see, e.g., [7] for details on the 
available algorithms and a recent review of numerical quadrature) are tailored 
for specific classes of integrands, with limited possibilities to solve simultane- 
ously families of integrals falling in different classes. 

We may therefore assume that a study able to increase the reliability of the 
automatic adaptive quadrature algorithms for solving parametric integrals in 
connection with the exploration of physical models is of interest for a great 
many users. Within the automatic adaptive quadrature, the approximate value 
Q of a given integral as well as its associated error estimate E are obtained as 
sums of local couples {q, e} of estimates over subranges. 

The general picture offered by the numerical evidence on the solution of para- 
metric integrals points towards the existence of a limited range of parameter 
values where the local quadrature sum q provides accurate solution of the in- 
tegral of interest, whereas for other parameter values the quadrature sum q 
is inaccurate. Over the range of accurate q outputs, the existing quadrature 
error estimators provide outputs e which, in most cases, grossly overestimate 
the actual quadrature error, whence the need of supplementary range subdi- 
visions and over computing to meet the input precision requests. However, 
over the range of inaccurate q outputs, the heuristics implemented in the local 
error estimators may result in spurious outputs quoted as reliable, hence the 
impossibility to detect such cases by means of the existing library codes. 

In the present paper we discuss a generalization of the approach proposed 
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in [8] intended to reconcile these two contradictory aspects. The cornerstone 
of such an analysis is the derivation of reliability criteria for the validation of 
the local error estimate e associated to a local quadrature sum q based on the 
study of the profile of the integrand at the set of quadrature knots entering 
the expression of q. 

The basic idea is that an unreliable estimate of e might originate either in the 
insufficient resolution of the integrand profile, or in the presence of difficult 
isolated points (integrable singularities, turning points, jumps) which result in 
slow convergence. The occurrence of each kind of difficulty can be evidenced 
by means of specific consistency criteria asking for the fulfillment of require- 
ments following from quite general considerations: the very definition of the 
Riemann integral, the fundamental properties of the basis polynomials which 
span the approximating linear space where the interpolatory polynomial of 
the quadrature rule is defined, the properties of the continuous functions at or 
near their extremal points, and the smoothness properties of the continuous 
functions inside their monotonicity subranges. 

If the integrand is well-conditioned but its profile is insufficiently resolved at 
the current set of quadrature knots, repeated subdivision of the integration 
range eventually results in the fulfillment of all the reliability constraints. A 
genuine difficult integrand point, however, recurs under repeated subrange 
subdivisions. Therefore, repeated analysis of the integrand profile under sub- 
range subdivision ultimately results in the diagnostics stability under iteration. 
This is the point where the general control routine of an automatic quadrature 
rule can take safe decisions concerning the best way to continue the solution 
refinement or to decide that the integral was solved within the input accuracy 
specifications. 

The paper starts with basic definitions and notations (section 2). In section 3, 
the main features of the validation procedure of a computed local couple {q, e} 
are discussed. Criteria for the identification of ill-conditioning features within 
an integrand profile are summarized in section 4. Their practical importance 
is assessed in the section 5 based on numerical evidence obtained from the 
solution of case study integrals by Gauss-Kronrod 10-21 quadrature rules [9] 
with improved error estimate [8]. Concluding comments are given in section 6. 



2 Definitions and notations 

Let / denote the actual value of the integral to be solved numerically, 

b 

I = I[f] = / g(x)f(x)dx , -oo < a < b < oo. (1) 
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Here, the weight function g(x) is an analytically integrable function which 
absorbs a difficult part of the integrand (e.g., an oscillatory or a singular 
factor). In the absence of such factors, g(x) = 1. The integrand function f(x) 
is assumed to be continuous almost everywhere on [a, b], such that (1) exists 
and is finite. 

A local quadrature rule produces as solution of (1) a couple {q, e}, where the 
quadrature sum q yields an approximate value of the integral I, while the local 
error estimate e > provides information on the accuracy of q. If e > j Cq | , 
where 

e Q = I-q (2) 

is the actual error associated to q, then the couple {q, e} is reliable, otherwise 
it is unreliable and the numerical solution fails. 

A (2n + l)-knot interpolatory quadrature sum q 2n is obtained as the analytical 
solution of the integral (1) with the integrand f(x) replaced by an interpola- 
tory polynomial of the 2ra-th degree, 

2n 

P2n(x) = a k Pk(x), (3) 
k=0 

where {^^(a;)} is the set of polynomials of degree at most In spanning the 
approximating space of P 2n (x). The coefficients are obtained from the set 
of conditions of interpolation 

P 2 n(x t ) = f( Xi ), (4) 

at a set of 2n + 1 abscissas (called quadrature knots) inside [a, b], 

a < x < xi < ■ ■ ■ < x 2n <b. (5) 

In the particular case of the symmetric (2n + l)-knot quadrature sums, the 
interpolation abscissas inside [a, b] are given by 

Xi = c + hz/i] c = (b + a)/2; h — (b — a)/2; i = —n, —n + 1, • • • , n, (6) 

where the reduced quadrature knots yi are defined on [—1, 1], such that = 
2/o < V\ < V2 < • • • < Vn < 1, while y_i = -y i: % = 1, • • • , n. 
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The local quadrature sum q2 n is then expressed as a linear combination of the 
integrand values at the quadrature knots, 



<?2n = <?2n[/] = W if( X i) > ( 7 ) 



with the quadrature weights showing the symmetry property W-i = Wi. 

The information provided by the In + 1 integrand values at the quadrature 
knots, {f(xi)\i = —n, ■ ■ ■ , n}, is insufficient for the derivation of an expression 
for the error estimate t2n associated to q2 n - 

Kronrod [10] derived an error estimate, called in what follows genuine Gauss- 
Kronrod (ggk) error estimate, from an upper bound of 



\Q2n-Qn\, (8) 



where q2 n is the quadrature sum (7), while q n is a lower degree quadrature 
sum derived over the subset of (6), 

n+7 ^ ra+7+2 < ' ' ' < £ n _^,_2 *C 7 • (9) 



Here, 7 = 1 for an open quadrature sum (typically, the Gauss-Kronrod (GK) 
quadrature where the spanning basis {pk(%)} in (4) is given by Legendre poly- 
nomials and their orthogonal Kronrod extensions), while 7 = for a closed 
quadrature sum (typically, the Clenshaw-Curtis (CC) quadrature where the 
spanning basis {pk(x)} i n (4) is given by Chebyshev polynomials). 

In what follows, the set of quadrature knots (6) is referred to as the fine dis- 
cretization of the integration domain [a, b], while the sparser set of quadrature 
knots (9) as the coarse discretization of [a, 6]. The integrand values at these 
knots define its fine and coarse samplings respectively. 

In the QUADPACK package [9], which has been incorporated in most major 
program libraries, while a ggk error estimate was implemented for the CC 
quadrature, the GK error estimate was reformulated as follows. Let / denote 
the computed value of the average of f(x) over [a, b] at the knots (6), 

f = q2n/(b-a), (10) 



and let A = Q 2n \f — f\ denote the computed value of J a fe \f(x) — f\dx, which 
measures the area covered by the deviations of f(x) around /. 
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The local QUADPACK error estimate (qdp) is then given by 
e qdp = A x min{(200e 9sfc /A) 3 / 2 , 1}. 



(11) 



The values (8) and (11) are taken for error estimates provided they do not fall 
below the attainable accuracy limit imposed by the relative machine precision. 
The latter threshold is defined as the product 



e ~roff = To€oQ 2 n \f\ 



;i2) 



Here r is an empirical multiplicative factor (following QUADPACK, we have 
chosen r = 50) and eo denotes the relative machine accuracy. 

For the case study integrals considered below, the value / of (1) is computed 
from the existing analytical expressions, such that the exact error e<g (2) of 
the quadrature sum q 2n can be defined. 

In the graphical representation of the quadrature errors, the moduli of the 
relative errors (simply called relative errors in the sequel) are useful, 

e a = \e a /I\, ae{2n,Q}. (13) 

The derivation of the local error estimates within a subroutine which imple- 
ments a quadrature rule uses information inferred from the estimated relative 
errors, 

Pa = \e a /q2n\ , ol G {ggk, qdp, 2n}. (14) 



3 Stability of the diagnostics under subrange subdivision 

Using (8), (11) and (12), we get the local error estimate [8] 

e 2n = max [e roff , mm(e ggk , e qdp )] , (15) 



the reliability of which is almost always subject to doubt, except for the case 
when the lower degree quadrature sum q n is sufficiently accurate such that 
the accuracy of q2 n itself reaches nine to ten significant figures at least. Such 
a condition can be confidently assumed to hold provided 

e 2 „ > e thr , e thr = 2~ 18 ~ 0.38 x lO" 5 . (16) 
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This empirically proposed threshold value is about two decimal figures more 
conservative than the smallest values of the unreliable computed error esti- 
mates over the evidence discussed in Sec. 5. 

If the opposite of (16) occurs, then a validation procedure is to be used to as- 
sess the reliability of the local couple {q2n, &2n}- Thus, a self-validating quadra- 
ture rule returns, besides the numerical output for e2 n , a flag having the value 
zero in case of assumed reliable outputs and non-zero value if the output is 
not validated. 

The validation procedure proposed in this paper is based on the study of the 
information contained in the profile of the integrand f(x) over [a, b], defined 
as the set of integrand values at the quadrature knots (5), completed with 
the endpoint values f(a) and f(b) in the case of open quadrature sums. Since 
the operation of subrange subdivision within automatic adaptive quadrature 
always involves inner abscissas at existing quadrature knots, the only price 
to be paid for the inclusion of the endpoint values in the integrand profile 
is the direct access of the general control routine to such data. This goal is 
achieved provided the generation of the integrand sampling at the quadrature 
knots (5) is done within a subroutine which is distinct from that implementing 
the quadrature rule and is directly subordinated to the general control routine. 

The study of the integrand profile starts with the definition of its monotonicity 
subranges, x^], over [a, b], where 

0, Xi < Xjj < X{ 2 < ■ • • < Xi m < Xi m+1 b , (17) 

denote the abscissas of the extremal points of f(x) within the sampling. 

In terms of the answer concerning the number of monotonicity subranges, 
several specific reliability criteria are checked and the number A of the in- 
fringements of these criteria is counted. There are three critical values of the 
pointer A in terms of which a decision is taken: 

• A = 0: probably the investigated integrand profile comes from a well- 
conditioned integrand, hence the output q2 n , Eq. (7) is reliable, while the 
quadrature error estimate (15) is overestimated. 

• A = lorA = 2: there is a high probability that a difficult isolated point is 
present which implies slow convergence of the quadrature sums. 

• A > 3: the insufficient resolution of the integrand profile at the involved 
quadrature knots is manifest. The output is useless and further subrange 
subdivisions are compulsory. 

The existence and finiteness of the Riemann integral (1) guarantees that, after 
a finite number of subrange subdivisions, the discretization process will reach 
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a stable profile configuration the refinement of which will result in unessential 
modifications only. 

Under the occurrence of isolated difficult points of the integrand, the dis- 
cretization process will resolve the profile over the well-conditioned subranges 
within a finite number of subrange subdivisions, and then it will mainly create 
a dense mesh around the difficult points. In this case, the automatic control 
subroutine will safely decide upon the activation of a specific convergence ac- 
celeration algorithm, such that a reliable numerical solution will be available 
in the end. 

The achievement of the stability of the diagnostics concerning the conditioning 
properties of the integrand profiles over subranges, got after a finite number of 
subrange subdivisions, is the fundamental feature which secures the efficiency 
of the procedure proposed in this investigation. 

The occurrence of consistent with each other reliability diagnostics over the 
current integration range and its subranges obtained by subrange subdivi- 
sion enables the general control routine to take safe decisions concerning the 
activation of the implemented alternative algorithms. 



4 Well-conditioned integrand profiles 

The consistency requirements satisfied by a well-conditioned integrand profile 
are formulated mostly locally and they follow from quite general considerations 
which are discussed in the next subsections. 

Any infringement of the consistency criteria derived below is to be added to 
the value of the ill-conditioning pointer A. 

4-1 Insensitivity of the integral sums to discretization details 

The standard definition of the integral sums in a Riemann integral assumes 
the fulfillment of the following two features: 

(i) The norm of the discretization step defined over the integration domain 
tends to zero. 

(ii) The integral sum is insensitive to the the addition or removal of a single 
discretization abscissa within the defined partition. 

In the quadrature algorithms, the norm of the discretization (6) is far from 
being close to zero. The quadrature knots are not distributed evenly either. 
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For the GK and CC quadrature rules mentioned above, the fundamental range 
[—1,1] consists of a sparser knot region centered around the origin and two 
denser knot regions located toward the range ends. The number of abscissas 
entering the integrand profile associated to a (2n + l)-knot open quadrature 
rule equals 2n + 3, while the corresponding number for a closed quadrature 
rule is 2n + 1. Therefore, for both kinds of quadrature rules, a particular 
inner reduced knot yi lies in the dense knot region provided the lengths of its 
two adjacent subranges are smaller than the threshold quantity for a uniform 
distribution, d av = 2/{2n + 3). 

An immediate consequence of the feature (i) is the property that the denser 
discretization regions of a smooth integrand f(x) secure better accuracy of 
their contributions to the quadrature sums than the sparser ones. We refor- 
mulate this observation as follows: the generation of the fine discretization (6) 
from the coarse discretization (9) is expected to result in non-essential modi- 
fications of the profile of f(x) over the regions of dense knot discretization. 

To characterize the extent to which a profile is modified by the addition of 
new knots inside the region of dense knot discretization, let us consider that 
Xq is such a knot. If x belongs to the set of extremal points (17) such that the 
integrand value f(xo) is isolated from the integrand values and f(x\) 

at the nearest neighbours x_i and x\ by the median line / = /, Eq. (10), then 
the knot xq is said to be sensitive. If both quantities f(x-i) and f(xi) stay 
on the same side with f(xo) with respect to the median line f — f, then the 
knot xq is said to be regular. If the median line f — f separates f(xo) from 
only one of the values f(x-i) or /(^i), then the knot x is said to be gray. 

We are now ready to formulate the first practical reliability criterion: 

(I) Non-sensitivity of the extremal points: 

The addition of supplementary quadrature knots to the coarse partition (9) 
to reach the fine partition (6) does not result in supplementary gray or 
sensitive extrema of the profile of a well-conditioned integrand over the 
regions of dense knot discretization. 

4-2 Features which stem from the basis polynomials 

Since the equations (6) perform the mapping of the original interval [a, b] 
onto the reduced interval [—1,1] over which the orthogonal polynomials are 
usually defined, in this subsection we refer to this reduced interval and use the 
notation pk(y) for the basis polynomials. All the properties discussed below 
hold true over arbitrary interval lengths, hence reference to the expression (3) 
of the interpolatory polynomial spanned by the basis orthogonal polynomials 
does not give rise to any confusion. 
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The existence and uniqueness of the interpolatory polynomial (3) is secured 
provided the set of basis polynomials spanning (3) define a Chebyshev system 
over [a,b]. Therefrom the following properties hold true: 

(iii) p (y) = const. 

(iv) The set of the successive extremal values of a polynomial Pk{y) of degree 
k > 1 defines an alternating sequence over [—1, 1]. 

(v) The zeros of the polynomials Pkiy) and p k+ i{y) are interlaced inside the 
open range (—1, 1). 

The average value /, Eq. (10), of the integrand f{x), which defines its zeroth 
order moment over the sampling (5) and is related to the coefficient of po(y) 
within a basis set of orthogonal polynomials, serves as reference value with 
respect to which the oscillations of the integrand profile are counted. The 
intersections of the integrand profile with the line f — f define the zeros of 
the integrand profile. 

The alternation property (iv) results in the important consequence that the 
deviations of the successive extremal values of a well-conditioned integrand 
profile from / define an alternating sequence with comparable amplitudes at 
the adjacent extremal knots (17). This property can be detailed for practical 
purposes into two well-conditioning alternation criteria: 

(Ha) Type-1 alternation criterion: 

- Each inner monotonicity subrange of a well- conditioned integrand profile 

intersects the line f — f. 

- The two end point monotonicity subranges do not diverge from f — f. 
(lib) Type-2 alternation criterion: 

Each inner gray extremal point which satisfies the type-1 alternation crite- 
rion is to stay sufficiently far from the line f = f. 

The test for the occurrence of an infringement of the type-1 alternation crite- 
rion is obvious. As it concerns the the latter criterion, two infringements are 
to be simultaneously tested: 

• The distance from f(x ) to / is to be smaller than those of its nearest 
neighbouring extrema. 

• Let a , di and a r denote the areas surrounded by / = / and the integrand 
profile around x Q and its nearest neighbours in the set (17). Then 

|a | < ti\ai + a r \, t 1 = 10, (18) 

where the value of t\ was chosen such as to point to a discrepancy exceeding 
an order of magnitude. The computation of the three local areas is done by 
compound trapeze rule which is robust and sufficiently accurate for the 
involved comparison. 
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Corroboration of the interlacing property (v) with the non-sensitivity criterion 
(I) results in a criterion for the distribution of the zeros of the integrand profile: 

(III) Non-sensitivity of the zeros: 

Over the dense knot regions, the numbers of zeros of the fine and coarse 
profiles of a well-conditioned integrand are the same. 



4-3 Integrand variations around its isolated extremal points 

The lateral first order derivatives of a smooth first order different iable function 
vanish at an extremal point, while the curvature of a second order differentiable 
function (which is given by the second order derivative) keeps constant sign 
over a nonvanishing neighbourhood of the extremum. 

Within the discrete mesh defined by the quadrature knots, inquiries about 
these properties can be made only at integrand profile approximations of iso- 
lated extremal points of the integrand. If xq is such a point, then a sufficiently 
large neighbourhood {£;,£ r } around xq can be defined within which the eval- 
uation of the quantities of interest is expected to be weakly influenced by the 
presence of neighbouring extrema. 

Let us assume that an isolated extremal point of a well-conditioned integrand 
was identified within a sufficiently well resolved integrand profile. The follow- 
ing consistency criteria establish well-conditioned behaviours of the data: 

(IV) First lateral derivative criterion: 

The approximation of the lateral first order derivatives at an isolated ex- 
tremum of the profile using hne sampling data is closer to zero as compared 
to the value estimated from data defined over a coarse sampling with respect 
to the extremum location. 
(V) Curvature sign constancy criterion: 

The sign of the second order derivative computed from fine sampling data 
centered at the extremum is the same as that of the value estimated from 
data involving the reference extremum as a lateral point to the left/right. 

We shall illustrate the quantitative implementation of these criteria for a ref- 
erence extremum xq which is said to be isolated to the right. That is, the 
neighbourhood {£z,£r} contains inside it the set of abscissas {x_±, Xo, xi, X2} 
at which the integrand function takes respectively the values /o, fi, fi\. 

To estimate the approximation of the first order right lateral derivative, we 
define the interpolatory polynomial of the third degree which fits these four 
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data. This yields the following result: 



tijine(xo) = d$ " ^ \h ,-ld% + . (19) 



h 2 ,-i 



Here, h^j = Xi—Xj, 42 = (fi—fj) / denote the first order divided differences 

at Xi and Xj, while 4^i = (42 ~~ 42) /^2,i and 4 2 -i = (42 — 4)2- 1) f^i-i 
denote specific second order divided differences. 

On the other hand, the coarse sampling around x yields: 

fr,coarse{x»)=d%. (20) 

The criterion (IV) then simply states that the approximations (19) and (20) 
should satisfy \f rJine {x G )\ < \ f' rtCoarae {xo)\. 

Over the same set of data, the criterion (V) requirement of constancy of the 
sign of the second order derivative results in the condition 

(42-42) (42-<y>o. (2i) 

For the extremal point x isolated to the left, similar conditions are derived 
from the data set {f-2, f-i, fo, fi} obtained at the abscissas {x_2, X-i,xq, x±}. 



4-4 Well- conditioning inside monotonicity subranges 



Inside any monotonicity subrange of a smooth first order differentiable func- 
tion f(x), the first order derivative varies smoothly from point to point. 

Within numerical quadrature, the fulfillment of this property for an integrand 
sampling can be checked by making use of first order divided differences. If 
the integrand profile is monotonic over [a, b], or monotonicity subranges can 
be defined which extend over three successive knots at least, then a smoothly 
varying profile will by characterized by the absence of jumps: 

(VI) Absence of jumps inside monotonicity subranges: 

Inside a monotonicity range, the ratio of two successive first order divided 
differences cannot exceed a relative smoothness threshold. 

If one of the knots involved in the divided differences is an extremal point, then 
this smoothness condition is to be checked only one-directionally, skipping the 
case of vanishingly small divided difference at the extremal point. 
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For knots far from inflection points, a threshold value tj mp = 10, correspond- 
ing to the agreement of the successive divided differences within an order of 
magnitude, is appropriate. In the neighbourhood of inflection points charac- 
terized by a maximum of the first order derivative, this value is to be halved 
to detect ill-conditioned behaviour, while in the neighbourhood of inflection 
points characterized by a minimum of the first order derivative, five times 
larger threshold value is appropriate. 



5 Numerical results 



The significance of the conditioning criteria discussed in the previous section 
is intuitive and straightforward. In addition to the case specified by the con- 
dition (16), a second case when the reliability analysis can be skipped is that 
of a monotonia profile characterized by an error estimate 

e 2n > 0.5|g 2n | . (22) 



Then the computed quadrature sum is highly inaccurate, such that an error 
flag can be directly assigned. 

The diagnostics of the reliability criteria (lib), (IV), (V), and (VI) depend on 
specific adjustable parameters. If the quantitative thresholds entering these 
criteria are decreased, the diagnostics becomes less permissive, with the con- 
sequence that the reliability range shrinks and the number of wrong diagnos- 
tics is decreased. The opposite occurs under the increase of the quantitative 
thresholds. The numerical data reported in this section show that, when cor- 
roborated with the requirement of the stability of the diagnostics formulated 
in Sec. 3, the formulation of the reliability criteria in Sec. 4 is able to elimi- 
nate practically all the spurious outputs occurring in an automatic adaptive 
quadrature algorithm. 

To illustrate the present analysis, a comparison is done of three codes using 
Gauss-Kronrod 10-21 (GK 10-21) quadrature rules: (a) the QUADPACK code [9], 
(P) the self- validating code of [8], and (7) the code using the present reliability 
analysis. 

Each code solved the parametric families of elementary integrals considered 
in ref. [8]. 
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The first is the family of integrals over [0, 1] of the terms of the fundamental 



power series, x , 



[ x n dx = — *— , n = 0, 1, • • • , 1023. 
J n + 1 



(23) 
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Fig. 1. Relative errors p q d p and p2 n , Eq. (14), eq, Eq. (13), of the GK 10-21 outputs 
for the family of integrals (23) at exponents n < 200. The upper leftmost solid 
line arrow points to the accuracy basin extension established by the QUADPACK code 
using the error estimate (11). The solid line arrow on the same vertical points to the 
upper accuracy of the quadrature sum q2 n retained as reliable by the QUADPACK code. 
The next pair of solid arrows show the result of the analysis done in ref . [8] . The left 
interrupted line arrow represents the extension of the reliability basin established 
by the present analysis, while the right one shows the exponent threshold above 
which the criterion (22) supersedes the need of reliability analysis. 

The integrands are monotonic, inflection points are absent over the integration 
range. Fig. 1 illustrates the behaviour of the error estimates with the power 
n running over the range {0,200}. The results obtained for this family of 
integrals can be summarized as follows: 

• The QUADPACK code infers an accuracy basin of the GK 10-21 code extending 
from n = to n max = 40, with the consequence that all the q^n outputs 
showing an actual accuracy lower than 14 decimal digits are thrown away. 
As shown in [8], this early cut of the accuracy basin does not rule out the 
possibility of wrong error diagnostics at asymptotically large n. 
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• The self- validating analysis of ref [8] extends the accuracy basin of GK 10-21 
up to n max = 59, which corresponds to a correct identification of the outputs 
q 2n as reliable up to accuracies of roughly nine significant digits. Above 
n > 60, all the reliability diagnostics have been correct. 

• The present analysis establishes an accuracy basin up to n max = 160, which 
corresponds to outputs q 2n showing at least three significant decimal digits. 
At exponents n > 195, the criterion (22) directly establishes the occurrence 
of unreliable q 2n outputs without making recourse to the reliability analysis. 

The second family solves integrals for a same integrand (which simulates a 
centrifugal potential at large x) over ranges of variable length, 

) 1 

/— dx = arctan(fr) , b = n, n = 0, 1, 10000. (24) 
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Fig. 2. Same as fig. 1 for the family of integrals (24), at upper integration ranges 
b < 260. The arrows bear the same significance. 

The integrands are monotonic, an inflection point is present. Fig. 2 illustrates 
the behaviour of the error estimates with the upper integration range b = n 
for n running over the range {0, 260}. In this figure the occurrence of cusps in 
the Eq curve points to the existence of fractional integration domain lengths 
at which the quadrature sum q 2n solves exactly the integral (24), such that 
the exact error changes sign. The results obtained for this family of integrals 
can be summarized as follows: 
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The QUADPACK code infers an accuracy basin of the GK 10-21 code extending 
up to n max = 10, with the consequence that all the q2 n outputs showing an 
actual accuracy lower than seven decimal digits are thrown away. All the 
QUADPACK reliability diagnostics above b = n = 2460 have been false. 
The self-validating analysis of ref [8] extends the accuracy basin of GK 10-21 
up to n max = 25, which corresponds to a correct identification of the outputs 
q2n as reliable up to accuracies of about six significant digits. At n > 26, all 
the reliability diagnostics have been correct. 

The present analysis establishes an accuracy basin up to n max = 37, which 
corresponds to outputs q 2n showing at least three significant decimal digits. 
At exponents n > 247, the criterion (22) directly establishes the occurrence 
of unreliable q2 n outputs without making recourse to the reliability analysis. 
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Fig. 3. Outputs of the GK 10-21 quadrature rule for the family of integrals (25) at 
p = 1. The significances of the solid line arrows are the same as in fig. 1. The left 
interrupted line arrow represents the extension of the reliability basin up to which 
the present analysis validates all the outputs q2 n - Inbetween the two interrupted 
line arrows, the diagnostics of the present analysis is too conservative in about one 
third of the solved cases. 

Next, we considered two pairs of families of integrals showing nonmonotonic 
(oscillatory) behaviour, written in algebraically equivalent forms: 



(CI) J e p(x - Xo) cos(cja:) dx = 



(25) 



16 
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Fig. 4. Same as fig. 3 for the family of integrals (28) at p = 1. 
i 

(C2) J 2e~ pXt> cosh(px) cos(cj;r) dx = (26) 
o 

= 2e' pxo [psinh(p) cos(cj) + cjcosh(p) sin(uj)}/(uj 2 + p 2 ) ; (27) 



(SI) J e p{x - Xa) sin(cux) dx = 



(S2) y 2e- p:co sinh(px) sin(c;x) 



= 2e px °[pcosh(p) sin(cj) — cusinh(p) cos(a;)]/(a; +p ) . 



(28) 

(29) 
(30) 



The parameter uj was chosen to run over the set of values 

u n = nn/60, ne {0,6000}, 



(31) 



while constant values p — 1 and x = — 1 have been chosen on the ground that 
they are typical for the description of the behaviour of the numerical results. 

The analysis of the families of integrals (25-29) shows that the identification of 
a well-conditioned nonmonotonic integrand profile needs testing the complete 
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Fig. 5. Same as fig. 3 for the family of integrals (26) at p = 1. 

set of consistency criteria established in Sec. 4. Therefore, the analysis is long. 
However, it is straightforward and can be easily implemented in a computer 
program. 

Figures 3 to 6 show outputs for the parameter n running over the range n G 
{0, 1080}. In these figures, two peculiarities of the Eq curves are apparent. 
Similar to Fig. 2, the occurrence of cusps at minima in the Eq curves point 
to the existence of values of the parameter u at which the given integrals are 
solved exactly by the quadrature sum q 2n , such that the exact error changes 
sign. The sharp maxima noticed in the Eq curves occur at u values which 
correspond to entire periods of the oscillatory factors over the integration 
range, such that important cancellation by subtraction effects occur which 
result in sensible worsening of the numerical output. 

A summary of the results obtained for the families of integrals (25), (26), (28), 
and (29) is given in Table I. 

The QUADPACK code predicts the narrowest accuracy basins in all the cases. 
Practically, any computed output qm with actual accuracy above the computer 
roundoff is ruled out as unreliable. At large values of the argument u of the 
trigonometric functions, this code results in an average rate of spurious outputs 
of about two percent. In figs. 3 and 4, unreliable estimates of this code are 
noticed at arguments u > 12n and u> > 147T respectively. The user is not 
notified of the wrong diagnostics associated to these outputs at the moment 
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Fig. 6. Same as fig. 3 for the family of integrals (29) at p = 1. 

of solving the integrals of interest. As mentioned in the Introduction, the only 
way of identifying them is the far end prediction of nonphysical results for the 
involved observables. 

Table 1 



Comparison of the stability basins and diagnostics reliability of the three codes 





Extension n max 
of the accuracy basins *) 


Number of spurious 
diagnostics at output **** 


Integral family 


(CI) (SI) (C2) (S2) 


(CI) (SI) (C2) (S2) 


QUADPACK 


219 192 359 379 


254 33 115 61 


Ref. [8] 


233 202 417 389 





Present 


399 395 820 769 
(521) (523) (1052) (1069) 


26 46 41 7 
(0) (0) (0) (0) 



*) For the present analysis, the upper values correspond to the left interrupted 
line arrows in Figs. 3 to 6. The values under parentheses correspond to the right 
interrupted line arrows in the same figures. 

**) For the present analysis, the upper values show the number of primary analysis 
failures. The vanishing values under parentheses show that all the primary analysis 
failures were corrected under subrange subdivision. 
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The self- validating procedure developed in ref. [8] slightly enlarged the exten- 
sion of the accuracy basin predictions, with no wrong outputs at all. 

The present reliability analysis identified substantially larger accuracy basins 
of the output. All the q 2n outputs showing more than six accurate figures have 
been correctly identified as reliable. For outputs q2 n showing inbetween six 
and three accurate figures, the present diagnostic was too conservative for 19 
(CI) integrals, 51 (SI), 57 (C2), and 176 (S2) integrals. At values of the 
argument u in large excess to those falling in the accuracy basins, a number 
of spurious diagnostics was produced by the primary reliability analysis. All 
the wrong diagnostics occurring at a first run were identified as wrong and 
corrected under subrange subdivision. 




-1 -0.5 0.5 1 -1 -0.5 0.5 




-1 -0.5 0.5 1 -1 -0.5 0.5 1 



Fig. 7. Ill-conditioned integrand features in the family of integrals (25-29) at 
u = 1612vr/60. 

Fig. 7 and fig. 8, show integrand profiles for the integrals (25-26) and (28- 
29) at the large parameter values oo = (16127r/60) and uo = (36467r/60) re- 
spectively, together with hints (showed by arrows) on infringements of the 
reliability criteria established in Sec. 4. 

A scrutiny of the integrand profiles shows that, in general, it is hardly probable 
that a highly oscillatory integrand structure can be resolved at the existing 
quadrature knots. However, if intermediate unresolved structures are present, 
these induce, as a rule, infringements of one or more reliability criteria. The 
complete list of criteria infringements is given below: 
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Fig. 8. Ill-conditioned integrand features in the family of integrals (25-29) at 
uj = 3646vr/60. 



• Criterion (I): 

In fig 7: the integral (CI) at and y 8 ; the integral (S2) at y 8 . 

In fig 8: the integral (SI) at y_io and y w (not shown in the plot); the integral 

(C2) at y 7 , and the integral (S2) at yio (not shown in the plot). 

• Criterion (11a) - over end subranges: 

In fig 8: the integral (SI) over the subrange [y_n,y_i ] (not shown in the 
plot). 

• Criterion (11a) - over inner subranges: 

In fig 7: the integral (CI) over the subranges [y_ 8 ,y_ 7 ], [y_ 7 ,y_ 6 ], [2/6,2/7], 
and [2/7,2/8] ; the integral (S2) over the subranges [y-7,y-e] and [y_ 6 ,y_ 5 ]; 
In fig 8: the integral (CI) over 15 (!) subranges: [y_i , 2/- 9 ], [y_ 9 ,y_ 8 ], 
[y-s,y-r}, [y-7,y-e\, [y_ 6 ,2/- 5 ], [y-4,y-s], [y- 3 ,Z/-2], [y-i,yo], [2/2,2/3]; [2/3,2/4], 
[2/5, 2/e], [2/6, 2/7], [2/7, 2/s], [2/8, 2/9], and [y 9 , y 10 ] (not shown in the plot since they 
are obvious). 

• Criterion (lib): 

In fig 7: the integral (S2) at y 7 ; 

In fig 8: the integral (SI) at y_ 9 and 2/9; the integral (C2) at y_6 and 2/3; 
the integral (S2) at y 9 . 

• Criterion (III): 

In fig 7: the integral (C2) inside the subranges: [y_g, y_ 8 ], and [?/_§, 2/— 7] ; 
the integral (S2) inside the subranges: [2/7,2/8], and [y 8 ,y 9 ]; 
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• Criterion (IV): 

In fig 8: the integral (SI) at y 9 (left derivative). 

• Criterion (V): 

In fig 8: the integral (S2) at y_ 6 (left neighbourhood). 

• Criterion (VI): 

In fig 7: the integral (SI) to the right of the knot y_2 and to the left of the 
knot 1/2; the integral (C2) to the right of the knot y_ 2 and to the left of the 
knot y 2 - 



6 Comments and conclusions 

The present investigation started from the need to get reliable numerical solu- 
tions of difficult parametric integrals occurring in theoretical models devoted 
to the study of the mechanism of the high-T c superconductivity in cuprates 
[l]-[5] and in a theoretical model of nuclear fission [6]. An important prereq- 
uisite to be satisfied by the automatic quadrature algorithms needed for the 
evaluation of the observables was the substantial increase of the reliability of 
the local error estimates. 

We have found that the study of the conditioning of the integrand profile 
enables the formulation of validation criteria (consistency conditions for a 
well-conditioned profile) able to identify insufficient profile resolution or the 
occurrence of isolated difficult points of the integrand. The analysis is simple, 
it is intuitive, it is easily implemented in a computer program and it is easily 
done. 

An important supplementary bonus offered by this analysis was the identifi- 
cation of q2n output reliability ranges which are substantially larger in com- 
parison with those obtained within the usual implementations of quadrature 
routines. The unsatisfactory features noticed in the validation criteria devel- 
oped in ref. [8] have been fully removed. 

The subroutines doing the profile analysis described in this paper are docu- 
mented and described in a separate document [11]. 

We conclude this study with the observation that the validation analysis de- 
scribed in the present paper is not intended to replace the existing quadrature 
algorithms. When the estimated accuracy exceeds a critical threshold (ten- 
tatively set to five decimal figures), then the present procedure is skipped 
altogether. However, if this threshold is not attained, it is automatically acti- 
vated by the general control routine. Its results prove to be invaluable in the 
analysis of complex integrands, where it is able to discover the overwhelming 
fraction of peculiar integrand profiles at early stages of the analysis. 
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